↓↓ Accès au projet ↓↓
Lien GitHub du projet


# Introduction


Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des operation statistique avancée grâce au langage R, à partir d’un jeu de données. Ainsi, notre choix s’est porté sur des données statistiques relatives à je sais pas. Celles-ci proviennent du site internet je sais pas.

Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer des régression linéaire sur notre jeu de données et produire un rapport.


# Outils et environnement de travail

R et RStudio

R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.

La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.


R markdown

L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.


Outils complémentaires




Données

Chargement des jeux de données :

Notre jeu de données comporte 102 individus décrits par 6 variables :

data <- read_csv("data/prestige.csv") %>% as_tibble()
## Parsed with column specification:
## cols(
##   education = col_double(),
##   income = col_double(),
##   women = col_double(),
##   prestige = col_double(),
##   census = col_double(),
##   type = col_character()
## )
data
## # A tibble: 102 x 6
##    education income women prestige census type 
##        <dbl>  <dbl> <dbl>    <dbl>  <dbl> <chr>
##  1      13.1  12351 11.2      68.8   1113 prof 
##  2      12.3  25879  4.02     69.1   1130 prof 
##  3      12.8   9271 15.7      63.4   1171 prof 
##  4      11.4   8865  9.11     56.8   1175 prof 
##  5      14.6   8403 11.7      73.5   2111 prof 
##  6      15.6  11030  5.13     77.6   2113 prof 
##  7      15.1   8258 25.6      72.6   2133 prof 
##  8      15.4  14163  2.69     78.1   2141 prof 
##  9      14.5  11377  1.03     73.1   2143 prof 
## 10      14.6  11023  0.94     68.8   2153 prof 
## # ... with 92 more rows


viz_edu <- data %>% 
  ggplot(aes(x = prestige, y = education, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(6, 16, 2)
  )

viz_edu <- ggplotly(viz_edu)

viz_inc <- data %>% 
  ggplot(aes(x = prestige, y = income, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(0, 20000, 10000)
  )

viz_inc <- ggplotly(viz_inc)

viz_women <- data %>% 
  ggplot(aes(x = prestige, y = women, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
      breaks = seq(25, 75, 25)
    ) +
    scale_y_continuous(
      breaks = seq(0, 100, 25)
    )

viz_women <- ggplotly(viz_women)

viz_census <- data %>% 
  ggplot(aes(x = prestige, y = census, col = type)) + 
  geom_point() + theme_bw() +
  scale_x_continuous(
    breaks = seq(25, 75, 25)
  ) +
  scale_y_continuous(
    breaks = seq(2500, 7500, 2500)
  )
  
viz_census <- ggplotly(viz_census)

subplot(viz_edu, viz_inc, viz_women, viz_census, nrows = 2, margin = 0.07)

Régression linéaire

En statistiques et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.

Dans un premier temps, nous pouvons étudier rapidement d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.

generate_scatterM <- function(data, mapping){
  data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) + 
    geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
    geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}

scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))

scatterM <- ggplotly(scatterM)
scatterM

Avec les observations précédentes, nous pouvons déjà écarter certaines linéarités et choisir lesquelles sont intéressantes à étudier. Ici, nous reprenons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable education (qui reflète le niveau d’étude). Cette fois nous utilisons ggplot pour une meilleure représentation.

education <- data %>% 
  ggplot(aes(x = education, y = prestige)) +
  geom_point() +
  geom_smooth(method = loess, se = T) +
  geom_smooth(method = lm, color = "red", se = F) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(6, 16, 2)
  ) +
  scale_y_continuous(
    breaks = seq(20, 80, 20)
  )

scatter_edu <- ggplotly(education)
scatter_edu

La droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.

La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.

On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.


Nous testons aussi la relation entre les variables income et prestige.

income <- data %>% 
  ggplot(aes(x = income, y = prestige)) +
  geom_point() +
  geom_smooth(method = loess, se = T) +
  geom_smooth(method = lm, color = "red", se = F) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(0, 25000, 5000)
  ) +
  scale_y_continuous(
    breaks = seq(20, 100, 20)
  )

scatter_inc <- ggplotly(income)
scatter_inc

Ici, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.

Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”.

lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu
## 
## Call:
## lm(formula = prestige ~ education, data = data)
## 
## Coefficients:
## (Intercept)    education  
##     -10.732        5.361

« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.

L’équation, de notre droite est donc \(y = 5.361x + -10.732\).

bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))

lag_plot <- bacfdf %>%
  ggplot(aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(aes(xend = lag, yend = 0)) +
  geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +     
  geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
  theme_bw()

lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot
durbinWatsonTest(lin_model_edu)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.2752512       1.43917       0
##  Alternative hypothesis: rho != 0
qq_edu <- lin_model_edu %>%
  ggplot(aes(sample = rstandard(.))) +
  stat_qq_line(color = "red", linetype = "dashed") +
  stat_qq(size = 1.2) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(-2, 2, 1)
  ) +
  scale_y_continuous(
    breaks = seq(-3, 3, 1)
  ) +
  labs(title = "Normal Q-Q",
       x = "Theoretical Quantiles\nlm(prestige ~ education)", 
       y = "Standardized residuals")

qq_edu <- ggplotly(qq_edu)
qq_edu
shapiro.test(residuals(lin_model_edu))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406
lin_model_inc <- lm(prestige~income, data = data)
lin_model_inc
## 
## Call:
## lm(formula = prestige ~ income, data = data)
## 
## Coefficients:
## (Intercept)       income  
##   27.141176     0.002897
qq_inc <- lin_model_inc %>%
  ggplot(aes(sample = rstandard(.)))+
  stat_qq_line(color = "red", linetype = "dashed") +
  stat_qq(size = 1.2) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(-2, 2, 1)
  ) +
  scale_y_continuous(
    breaks = seq(-3, 3, 1)
  ) +
  labs(title = "Normal Q-Q",
       x = "Theoretical Quantiles\nlm(prestige ~ income)", 
       y = "Standardized residuals")

qq_inc <- ggplotly(qq_inc)
qq_inc
shapiro.test(residuals(lin_model_inc))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lin_model_inc)
## W = 0.97169, p-value = 0.02729
sc_loc_edu <- lin_model_edu %>%
  ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) + 
  geom_point(size = 1.2) +
  geom_smooth(method = loess, se = FALSE) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(30, 70, 10)
  ) +
  scale_y_continuous(
    limits = c(0, 1.5),
    breaks = seq(0, 1.5, 0.5)
  ) +
  labs(title = "Scale-location",
       x = "Fitted values\nlm(prestige ~ education)", 
       y = "sqrt |Standardized residuals|")

sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_edu
ncvTest(lin_model_edu)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.6327545, Df = 1, p = 0.42635
res_v_fit_edu <- lin_model_edu %>%
  ggplot(aes(fitted(.), residuals(.))) + 
  geom_point() +
  stat_smooth(method="loess", se = FALSE) + 
  geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(30, 70, 10)
  ) +
  scale_y_continuous(
    limits = c(-30, 20),
    breaks = seq(-30, 20, 10)
  ) +
  labs(title = "Residuals vs. Fitted values",
       x = "Fitted values\nlm(prestige ~ education)", 
       y = "Residuals")
    
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_edu
summary(lin_model_edu)
## 
## Call:
## lm(formula = prestige ~ education, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0397  -6.5228   0.6611   6.7430  18.1636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -10.732      3.677  -2.919  0.00434 ** 
## education      5.361      0.332  16.148  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared:  0.7228, Adjusted R-squared:   0.72 
## F-statistic: 260.8 on 1 and 100 DF,  p-value: < 2.2e-16
confint(lin_model_edu)
##                  2.5 %    97.5 %
## (Intercept) -18.027220 -3.436744
## education     4.702223  6.019533
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence")
##        fit      lwr      upr
## 1 38.64170 36.58966 40.69374
## 2 48.61293 46.81134 50.41452
## 3 67.64405 64.52387 70.76423
data$residuals <- residuals(lin_model_edu)
head(data)
## # A tibble: 6 x 7
##   education income women prestige census type  residuals
##       <dbl>  <dbl> <dbl>    <dbl>  <dbl> <chr>     <dbl>
## 1      13.1  12351 11.2      68.8   1113 prof       9.25
## 2      12.3  25879  4.02     69.1   1130 prof      14.1 
## 3      12.8   9271 15.7      63.4   1171 prof       5.67
## 4      11.4   8865  9.11     56.8   1175 prof       6.31
## 5      14.6   8403 11.7      73.5   2111 prof       5.86
## 6      15.6  11030  5.13     77.6   2113 prof       4.49
data$fitted <- fitted(lin_model_edu)
head(data)
## # A tibble: 6 x 8
##   education income women prestige census type  residuals fitted
##       <dbl>  <dbl> <dbl>    <dbl>  <dbl> <chr>     <dbl>  <dbl>
## 1      13.1  12351 11.2      68.8   1113 prof       9.25   59.5
## 2      12.3  25879  4.02     69.1   1130 prof      14.1    55.0
## 3      12.8   9271 15.7      63.4   1171 prof       5.67   57.7
## 4      11.4   8865  9.11     56.8   1175 prof       6.31   50.5
## 5      14.6   8403 11.7      73.5   2111 prof       5.86   67.6
## 6      15.6  11030  5.13     77.6   2113 prof       4.49   73.1
vcov(lin_model_edu)
##             (Intercept)  education
## (Intercept)   13.520978 -1.1835055
## education     -1.183505  0.1102162
pred_interval <- predict(lin_model_edu, interval = "prediction")
## Warning in predict.lm(lin_model_edu, interval = "prediction"): predictions on current data refer to _future_ responses
data <- cbind(data, pred_interval)
linear_regression <- data %>% 
  ggplot(aes(y = prestige, x = education)) +
  geom_point(size = 1.2) +
  geom_smooth(color = "red", method = "lm", fill = "red") +
  geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  geom_line(aes(y =  upr), color = "red", linetype = "dashed") +    
  annotate("text", x = 9, y = 80, 
           label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3), 
                          " + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
  theme_bw() +
  scale_x_continuous(
    breaks = seq(6, 16, 2)
  ) +
  scale_y_continuous(
    breaks = seq(25, 75, 25)
  ) +
  labs(title = "Linear Regression",
       x = "Education", y = "Prestige")
  

linear_regression <- ggplotly(linear_regression)
linear_regression